Statistically Steady Turbulence in Soap Films: Direct Numerical Simulations with 

Ekman Friction 
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We present a detailed direct numerical simulation (DNS) designed to investigate the combined 
effects of walls and Ekman friction on turbulence in forced soap films. We concentrate on the 
forward-cascade regime and show how to extract the isotropic parts of velocity and vorticity structure 
functions and thence the ratios of multiscaling exponents. We find that velocity structure functions 
display simple scaling whereas their vorticity counterparts show multiscaling; and the probability 
distribution function of the Weiss parameter A, which distinguishes between regions with centers 
and saddles, is in quantitative agreement with experiments. 
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The pioneering work of Kraichnan jlj showed that fluid 
turbulence in two dimensions (2D) is qualitatively differ- 
ent from that in three dimensions (3D): in the former we 
have an infinity of extra conserved quantities, in the in- 
viscid, unforced case; the first of these is the enstrophy. 
It turns out, therefore, that 2D turbulence displays an 
inverse cascade of energy, from the length scale at which 
the force acts to larger length scales, and a forward cas- 
cade of enstrophy, from the forcing length scale to smaller 
ones; by contrast, 3D turbulence is characterised by a for- 
ward cascade of energy [2[. Kraichnan's predictions were 
first confirmed in atmospheric experiments in quasi-two- 
dimensional, stratified flows Q; subsequent experiments 
have studied systems ranging from large-scale geophysi- 
cal flows to soap films [MH, H, @, 0, 01 ■ The latter have 
proved to be especially useful in characterizing 2D tur- 
bulence. 

We present the first direct numerical study (DNS) that 
has been designed specifically to explore the combined ef- 
fects of the air-drag induced Ekman friction a and walls 
on the forward cascade in 2D turbulence; and we employ 
the Kolmogorov forcing used in many soap-film exper- 
iments [H, 0, 0, S|- Thus we can make a far more de- 
tailed comparison with these experiments than has been 
attempted hitherto. In particular, if we use values of 
a that are comparable to those in experiments, we find 
that the energy dissipation rate because of the Ekman 
friction is comparable to the energy dissipation rate that 
arises from the conventional viscosity. We show how 
to extract the isotropic parts [§| of velocity and vortic- 
ity structure functions and then by using the extended 
self-similarity (ESS) procedure [10[, we obtain ratios of 
multiscaling exponents whence we conclude that veloc- 
ity structure functions show simple scaling whereas their 
vorticity counterparts display multiscaling. Most impor- 
tant, our probability distribution function (PDF) of the 
Weiss parameter A is in quantitative agreement with 
that found in experiments 0, 01 • 

For the low-Mach-number flows we consider, we can 



use the following soap-film equations 0, [l^ : 

{dt + n-V)uj ^ vV'^Lo - auj + F^l p; V'^tp^Lu. (1) 



Here 



{-dy%l},dxil}), ip, and lu 



V X u are, re- 



spectively, the velocity, stream function, and vorticity 
at the position x and time t; we choose the uniform 
density p ~ 1; a is the Ekman friction coefficient, h' 
is the kinematic viscosity, and F^^ = kinjFoCos{kinjx), 
a Kolmogorov- type forcing term, with amplitude Fq, 
and injection wave vector kinj (the length scale £inj = 
'^T^/kinj)- We impose no-slip {ip = 0) and no-penetration 
{Vip-n = 0) boundary conditions on the walls, where n is 
the outward normal to the wall. If we non-dimensionalize 



X by 
with 



t by k-, 



inr - ■"•mj/^^. and F^ by 27r/(fci„j-||F^||2), 
F^||2 = (£4 l-FLpdx)^/^, then we have two con- 
trol parameters, namely, the Grashof [l3| number Q = 
27i'||-FL;||2/(fcf„jPJ^^) and the non-dimensionalized Ekman 
friction 7 = 0;/ {kf^ji/). For a given set of values of G and 
7, the system attains a noncquilibrium statistical steady 
state after a time t/r ~ 2.8, where r L/urms is the 
box-size time, L the side of our square simulation do- 
main, and Urms the root-mean-square velocity. In this 
state the Reynolds number Re = Urms/{kinji^), the en- 
ergy, etc., fluctuate; their mean values, along with one- 
standard-deviation error bars, are given in Table U that 
lists the values of the parameters in our runs Rl-7. 

We use a fourth-order Runge-Kutta scheme with step 
size 5t = 10~^ for time marching in Eq. ([1]) and evalu- 
ate spatial derivatives via second-order and fourth-order, 
centered, flnite differences, respectively, for points adja- 
cent to the walls and for points inside the domain. The 
Poisson equation in ^ is solved by using a fast-Poisson 
solver 1^ and uj is calculated at the boundaries by using 
Thom's formula [l^. To evaluate spatiotemporal aver- 
ages, we store '0(x, t„) and uj{x, tn), with i„ = (4-|-nA)r, 
n 0, 1, 2, . . . , n^ax: and 96 < rimax < 200; A = 0.28 
for runs Rl-6 and A = 0.13 for run R7. 

Figures [T](a)-(f) show the time evolution of the kinetic 
energy E{t) = (J^u^dx)/^, viscous energy-dissipation 
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Rl 


1025 


0.45 


45 


0.25 


3.5 


23.3 ± 0.4 


15.1 ± 0.5 


-28 ± 2 


-13.6 ± 0.5 


5.3 ± 0.3 


0.32 ± 0.01 


3.1 ± 0.1 


R2 


1025 


1.25 


45 


0.71 


3.5 


19.6 ± 0.3 


10.7 ± 0.3 


-28 ± 1 


-26.8 ± 0.9 


4.8 ± 0.2 


0.33 ± 0.01 


3.1 ± 0.1 


R3 


1025 


1.25 


60 


0.71 


4.7 


24.0 ± 0.5 


15.9 ± 0.6 


-40 ± 2 


-39.9 ± 1.4 


7.2 ± 0.4 


0.33 ± 0.01 


2.8 ± 0.1 


R4 


2049 


0.45 


45 


0.25 


3.5 


23.2 ± 0.4 


15.1 ± 0.5 


-28 ± 2 


-13.6 ± 0.4 


5.3 ± 0.3 


0.31 ± 0.01 


3.2 ± 0.2 


R5 


2049 


1.25 


45 


0.71 


3.5 


19.6 ± 0.4 


10.8 ± 0.4 


-28 ± 1 


-27.0 ± 1.0 


4.8 ± 0.2 


0.33 ± 0.01 


3.1 ± 0.1 


R6 


2049 


1.25 


60 


0.71 


4.7 


23.8 ± 0.4 


15.9 ± 0.6 


-40 ± 2 


-40.0 ± 1.5 


7.2 ± 0.4 


0.33 ± 0.01 


3.1 ± 0.1 


R7 


3073 


0.45 


45 


0.25 


3.0 


26.5 ± 0.4 


20.0 ± 1.0 


-26 ± 2 


-17.8 ± 0.6 


5.0 ± 0.4 


0.31 ± 0.01 


3.7 ± 0.3 



TABLE I: Parameters for our runs Rl-7: N, the number of grid points along each direction, 7, Q, Re (we use 1/ = 0.016, linj = 
0.6, and a square simulation domain with side L — 7, grid spacing 5^ = L/N, area A, and boundary dA), the time-averaged 
kinetic energy, viscous-energy-dissipation rate, and the energy-dissipation rate because of Ekman friction, E, e^, and e^, respec- 
tively, A = y^iid^u'yy) {{dyU'^y) , b = -{d^u'ydyu'^}/A, and the boundary-layer thickness 5^ oc ((/^^ u'^/ /g^((Vtj).n)^)^''^> 
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FIG. 1: (Color online) Representatitve plots from runs Rl(red 
circles), R2(black lines), and R3(black squares), showing the 
time evolution of E{t)/ME [{a) and (6)], e^{t)/Af [(c) and 
(d)], and te{t)/M [(e) and (/)]. In (a), (c), and (e) we keep 
Q fixed and vary 7 (7 = 0.25 (red circles) and 7 = 0.71 (black 
line)). In (6), (d), and (/) we maintain Re ~ 21.2 by varying 
7 (7 = 0.25(red circles) and 7 = 0.71(black squares)) and Q. 



rate ev{t) 



/(J^ |ct;|^c?x)/A, and energy-dissipation 



rate because of the Ekman friction ee(t) 



-2aE(t) 



^4 ^.^3) 

= {ey{t)), and e 



and Af ^ -kf,^^-^ 



(non-dimensionalized, respectively, by TVb = {vkinjY 
The mean values E = {E{t)), 
ee{t)), given in Table |T1 are com- 
parable to those in experiments; note that ei, and eg are 
of similar magnitudes. By comparing data from runs Rl 
(red circles) and R2 (black hues) in Figs.[Ha), (c), and 
(e) we see that, if we fix Q and increase 7, E decreases, e^, 
remains unchanged (within error bars), and eg increases. 
If we change both Q and 7, we can keep the mean Re 
fixed, as in runs Rl and R3 in Table HI by compensat- 
ing an increase in 7 with an increase in Q (cf. Ref. [^); 
in Figs. [Ijb), (d), and (f) we see, by comparing runs 
Rl (red circles) and R3 (black squares), that E remains 
unchanged (within error bars), whereas both e^, and Cg 
increase as 7 and Q increase in such a way that Re is held 
fixed. 

Since Kolmogorov forcing is inhomogeneous, we use 
the decomposition ip ~ (tp) + tp' and cu = {us) + 



Lu' , where the angular brackets denote a time aver- 
age and the prime the fluctuating part [lOj, to calcu- 
late the order-p velocity and vorticity structure func- 
tions S'p(rc,R) = (|(u'(rc + R) -u'(r)) • R/i?|P) and 
^^(rcR) = (|t^'(rc + R) - t^'(r)|^'), respectively, where 
R has magnitude R and rc is an origin. Figures [2fa) and 
(b) show pseudocolor plots of S'2(rc,R) and S'2(rc,R), 
respectively, for Tc = (2, 2); other values of Tc yield sim- 
ilar results so long as they do not lie near the bound- 
ary layer (Table |T| of thickness Sb {vc is chosen at 
least bSh away from all boundaries). We now calcu- 
late 52(R) = (52(re,R))re and S^R) = (52"(re, R))r., 
where the subscript rc denotes an average over the origin 
(we use Tc = {i,j),2 < i,j < 5); these averaged struc- 
ture functions [Figs. [DJc) and (d)] are nearly isotropic 
for R < linj but not so for R > imj. To obtain the 
isotropic parts in an SO{2) decomposition of these struc- 
ture functions we integrate over the angle 9 that R 
makes with the x axis to obtain Sp{R) = 5p(R)d6' 

and S^{R) = J^'' S!^{Ii)de. Given Sp{R) and S^im 
we use the extended-self-similarity (ESS) procedure [10| 
to extract the multiscaling-exponent ratios Cp/C^ and 
Cp/C2i respectively, from the slopes (in the forward- 
cascade inertial range) of log- log plots of Sp{r) versus 
52 (r) [Fig. He)] and S^{R) versus S^{R) [Fig.^f)] [H. 
The insets Figs. [2] (el) and (fl) show, respectively, plots 
of the local slopes Xp = ^logio Sp{R)/d\ogiQ S2{R) ver- 
sus logio52(i?) and = d\og^^ {R) / d\og^^ {R) 
versus \ogiQS2{R) in the forward-cascade regime; the 
mean values of Xp and Xpi over the ranges shown, yield 
the exponent ratios Cp/G and C,p/C,2 that are plotted 
versus p in Figs. O (e2) and (f2), respectively, in which 
the error bars indicate the maximum deviations of Xp 
and Xp fi'om their mean values. The Kraichnan-Leith- 
Batchelor (KLB) predictions [l| for these exponent ratios, 
namely, C^^^/Cz"^^^ ^ r^'^ and (^-^^klb /(-^.klb ^ 
agree with our values for Cp/C2 but not Cp IO2' velocity 
structure functions do not display multiscaling [Fig. [2] 
(e2)] whereas their vorticity analogs do [note the curva- 
ture of the plot in Fig. [5] (f2)]. This is in consonance 
with the results of DNS studies with periodic boundary 



3 



conditions [13, 1S|- Indeed, if we use the same values 
of 7 as in Ref. 17[ , we obtain the same exponent ratios 
(within error bars); thus our method for the extraction of 
the isotropic parts of the structure functions suppresses 
boundary and anisotropy effects efficiently. 

For an inviscid, incompressible 2D fluid the local flow 
topology can be characterized via the Weiss criterion [llj 
that uses the invariant A = (w^ — (T^)/8, where = 



and 



{diUj + djUi 



This criterion pro- 



vides a useful measure of flow properties even if > as 
noted in the experiments of Ref. [6| : Regions with A > 
and A < correspond to centers and saddles as wc show 
in Fig. [3] (a) by superimposing, at a representative time, 
a pseudocolor plot of A on contours of V'- This result is 
in quali_t_ativc accord with experiments [see, e.g., Fig. 1 



of Ref. [6[ and also earlier DNS studies [ll| , which do not 
use Ekmaii friction]. In Fig.[3](b) we compare the scaled 
PDFs P2{A/ Arms) with data obtained from points near 
the walls (black curve) and from points in the bulk (red 
curve); the clear difference between these, not highlighted 
before, indicates that regions of large A are suppressed 
in the boundary layers. There is a generation of strain 
and vorticity in these boundary layers and scatter plots, 
not shown, indicate ^ here; this leads to the sup- 
pression of regions of large A in P2{A/ Arms)- 

Figures [3] (c) and (d) show the PDF Pi (A) and the 
scaled PDF P2{A/Arms) for runs R4 (red hue) and R5 
(blue dashed line), with 7 = 0.25 and 7 = 0.71, respec- 
tively, and Q = 3.5 x 10^; by comparing these figures we 
see that both Pi and P2 overlap within error bars for 
runs R4 and R5. We believe this is because, in fixed-^ 
runs like R4 and R5, e„ does not change [Table [l] even 
though 7 changes. By contrast, if we compare Pi and P2 
[Figs. [3] (c) and (f)] for runs R4 (red line) and R6 (blue 
dashed line) , in which the mean Re is held fixed by tuning 
both 7 and Q, wc find, in agreement with experiments 0, 
that the PDFs Pi do not agree for these runs, but the 
PDFs P2 overlap within error bars. Our results for P2 in 
Fig-[S{f) are in quantitative agreement with experiments: 
we have obtained the points in this plot by digitising data 



points [see http:/ /_www.fran tz.fi/software/g3data.php| in 
Fig. 2(d) of Ref. Q ; the errors in these points are com- 
parable to the spread of data in Q . Conditional expec- 
tation values of (cr^} and (w^), for a given value of A, 
also agree well with experiments as can be seen by com- 
paring Fig. m with Fig. 3 of Ref. 0. We also present 
in Figs. [4] (b-d) pseudocolor plots of the joint PDFs of 
duj{r) EE uj'{:x.+rex)-uj'{x), Su^ir) = u'^ix+re^) -u'^{x) 



or 



Suip{r) 



u'y{-x.) with A' = det(M), 
1^(3 ^ dau'g, H. a circular 



Uy{x. + re^) 
M^P ^ /^m"/5dr/^dr, 
disc with center at x -|- {r/2)ex and radius r/2, and r in 
the forward-cascade regime; we obtain striking agreement 
with experiments as can be seen by comparing Figs. [4] 
(c-d) with Figs. 1-2 of Ref. 0. Finally, we calculate 

'{i9.u'y)^){idyK)^) and b ^ - {d^u'^dyiQ / A (see 



TableHl) and obt ain excellent agreement with experiments 

Some earlier numerical studies of 2D, wall-bounded, 
statistically steady turbulent flows [l3| use forcing func- 
tions that are not of the Kolmogorov type; furthermore, 
they do not include Ekman friction. Other numerical 
studies, which include Ekman friction and Kolmogorov 
forcing, employ periodic boundary conditions [l7, 18, 1^. 
To the best of our knowledge our study of 2D turbulent 
flows is the first one that accounts for Ekman friction, 
realistic boundary conditions, and Kolmogorov forcing. 
Thus we can make quantitative comparisons with soap- 
film experiments; and the agreement between our results 
and those of Refs. 0, 0, 0, Q vindicates the use Eq. ([1]) 
as a model for these soap films We hope our re- 

sults will stimulate experimental studies designed to ex- 
tract (a) the isotropic parts of structure functions (and 
thereby to probe the multiscaling of vorticity structure 
functions [Fig.[2](fl)]) or (b) the PDF P2(A/A„„s) [Fig.[3] 
(b)] near soap-film boundaries. 
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computational facilities. RP is a member of the Interna- 
tional Collaboration for Turbulence Research. 
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FIG. 2: (Color online) Pseudocolor plots of (a) S2{ra, R), for rc = (2, 2), (b) S2{R) (average of 52(rc, R) over rc), (c) S'^(rc, R), 
for Tc — (2, 2), and (d) S2iR) (average of 52 (rc, R) over rc). Log-log ESS plots of the isotropic parts of the order-p (e) velocity 
structure functions Sp{R) versus S2{R) and (f) the vorticity structure functions Sp{R) versus S2 (J?); p = 3 (blue line with 
circles), p — 4 (green line with triangles), p — 5 (red line with squares), and p = 6 (cyan line with stars); plots of the local 
slopes Xp and Xp (see text), in the forward-cascade inertial range: (el) Xp versus log^g S2{R) and (fl) Xp versus logj^Q S2{R)- 
Plots versus p of the exponent ratios (e2) Cp/Ca, along with the KLB prediction (red line), and (f2) Cp /C2' and error bars from 
the local slopes (see text). All plots are for run R7. 
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FIG. 3: (Color online) (a) Representative pseudocolor plot 
of A superimposed on a contour plot of the stream function 
^; (b) the PDF -Pi (A) obtained from points in the bulk Si, < 
x,y < L — Sb (red line) and from those points within a distance 
5b from the boundaries (black line) for our run R7; plots of (c) 
Pi (A) versus A and (d) P2(A/Arms) versus A/ Arms for fixed Q 
and 7 = 0.25(red line) and 7 = 0.71(blue dashed line) [runs R4 
and R5]; plots of (e) -Pi (A) versus A and (f) P2 (A/Arms) versus 
A/Arms [runs R4 and R6 with Re ~ 23.5] and 7 = 0.25(red 
line) and 7 = 0.71 (blue dashed line) and points (black dots) 
extracted from Fig. 2(d) of Ref. 0|. In (c)-(f) the fluctuating 
part of the velocity is used for A. One-standard-deviation 
error bars are indicated by the shaded regions. 
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0.12), A' /A'rms) for our run R7. The contours and the shading 
are for the logarithms of the joint PDFs. 
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